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DIFFERENTIAL ROTATION IN NEUTRON STARS: MAGNETIC BRAKING AND VISCOUS DAMPING 
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ABSTRACT 

Diffferentially rotating stars can support significantly more mass in equilibrium than nonrotating or uniformly 
rotating stars, according to general relativity. The remnant of a binary neutron star merger may give rise to such 
a "hypermassive" object. While such a star may be dynamically stable against gravitational collapse and bar 
formation, the radial stabilization due to differential rotation is likely to be temporary. Magnetic braking and 
viscosity combine to drive the star to uniform rotation, even if the seed magnetic field and the viscosity are small. 
This process inevitably leads to delayed collapse, which will be accompanied by a delayed gravitational wave 
burst and, possibly, a gamma-ray burst. We provide a simple, Newtonian, MHD calculation of the braking of 
differential rotation by magnetic fields and viscosity. The star is idealized as a differentially rotating, infinite 
cylinder consisting of a homogeneous, incompressible conducting gas. We solve analytically the simplest case in 
which the gas has no viscosity and the star resides in an exterior vacuum. We treat numerically cases in which the 
gas has internal viscosity and the star is embedded in an exterior, low-density, conducting medium. Our evolution 
calculations are presented to stimulate more realistic MHD simulations in full 3+1 general relativity. They serve 
to identify some of the key physical and numerical parameters, scaling behavior and competing timescales that 
characterize this important process. 



1. INTRODUCTION 

Neutron stars may form with appreciable differential rota- 
tion. The merger of a binary neutron star system, for exam- 
ple, will form a differentially rotating remnant, an outcome 
considerably enhanced by the liklihood that the neutron stars 
are nearly irrotational before coalescence (Bildsten and Cutler 
1992; Kochanek 1992; Lai, Rasio & Shapiro 1993). The forma- 
tion of a differentially rotating remnant following binary merger 
was first demonstrated in Newtonian hydrodynamic simula- 
tions (see e.g., Rasio and Shapiro 1992,1994,1999) and dra- 
matically confirmed in fully relativistic simulations (Shibata & 
Uryu 2000), where the remnant's core is found to rotate ap- 
preciably faster than its envelope. Core collapse in a super- 
nova may also result in a differentially rotating neutron star, 
even if the rotation of the progenitor at the onset of collapse is 
only moderately rapid and almost uniform (see, eg, Zwerger & 
Miiller 1997; Rampp, Miiller & Ruffert 1998, and references 
therein). Conservation of angular momentum during collapse 
causes the ratio (3 = T/\W\, where T is the rotational kinetic 
energy and W is the gravitational potential energy, to grow like 
ST 1 , where R is the equatorial radius. Thus j3 can increase by 
several orders of magnitude during the implosion. Since uni- 
formly rotating, compressible stars can only support very small 
values of (3 in equilibrium without shedding mass at the equator 
(see, e.g., Tassoul 1978, Shapiro & Teukolsky 1983), collapsed 
cores of rapidly rotating progenitors must acquire some degree 
of differential rotation at birth as they settle into equilibrium. 
Similar arguments hold for accretion-induced collapse of white 
dwarfs to neutron stars. 

Differentially rotating neutron stars can support significantly 
more rest mass than their nonrotating or uniformly rotating 
counterparts. Baumgarte, Shapiro & Shibata (2000; hereafter 
BSS), have recently constructed relativistic equilibrium mod- 
els of differentially rotating "hypermassive" polytropes. These 
configurations have rest masses that exceed both the maximum 



rest mass of nonrotating spherical stars (the TOV limit) and 
uniformly rotating stars at the mass-shedding limit ("supramas- 
sive stars"), all with the same polytropic index. BSS performed 
dynamical simulations in full general relativity to demonstrate 
that hypermassive stars can be constructed that are dynamically 
stable against radial collapse and nonradial bar formation (see 
also Shibata, Baumgarte & Shapiro 2000). The fully relativis- 
tic binary coalescence calculation of Shibata & Uryu (2000) 
reveals that such a dynamically stable hypermassive star can 
arise following binary merger. BSS have argued that the dy- 
namical stabilization of a hypermassive remnant by differential 
rotation may lead to delayed collapse and a delayed gravita- 
tional wave burst. The reason is that the stabilization due to 
differential rotation, although expected to last for many dynam- 
ical timescales (i.e. many milliseconds), will ultimately be de- 
stroyed by magnetic braking and/or viscosity. These processes 
drive the star to uniform rotation, which cannot support the ex- 
cess mass, and will lead to catastrophic collapse, possibly ac- 
companied by some mass loss. 

The purpose of this paper is to provide a time-dependent cal- 
culation of the damping of differential rotation by magnetic 
braking and viscosity in a rapidly spinning star. The goal is 
to construct an exact, albeit idealized, Newtonian model which 
readily reveals some of the important physical effects, together 
with their timescales and scaling behavior. By adopting a suf- 
ficiently simple, laminar model with a high degree of spatial 
symmetry, we are able to perform a tractable analysis of the 
evolution from first principles and without approximation. Our 
motivation for this simplified approach is to help identify the 
minimum ingredients for a more realistic, future numerical cal- 
culation in full general relativity. We wish to demonstrate ex- 
plicitly why it is essential to include magnetic fields in fully 
relativistic simulations of the late stages of coalescing binary 
neutron stars and other dynamical scenarios where differential 
rotation naturally arises and significantly influences the struc- 
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ture and stability of the final configuration. The results have 
important implications for the detection of gravitational wave 
and gamma-ray bursts (BSS; Spruit 1999a). 

Our discussion is by no means the first, (see, e.g., Spitzer 
1978 for discussion and references on magnetic braking of an- 
gular momentum in interstellar clouds; see Spruit (1999b) and 
references therein for a discussion of differential rotation and 
the evolution of initially weak magnetic fields in stellar interi- 
ors), and surely will not be the last on this nontrivial topic, es- 
pecially given the simplified nature of our model. However, the 
importance of magnetic braking has been pointed out only re- 
cently in the context of binary neutron star coalescence, differ- 
ential rotation of the remnant and the emission of gravitational 
waves following delayed collapse (BSS). Presumably the most 
effective way to stimulate numerical relativists to include mag- 
netic fields in 3 + 1 fully relativistic hydrodynamics codes now 
under development is to provide a concrete calculation demon- 
strating their physical importance, which we do below. 

In Section 2 we describe our basic model. In Section 3 we 
set out the fundamental equations and put them in a convenient 
nondimensional form. In Section 4 we summarize solutions for 
three distinct cases, depending on whether the stellar exterior 
region is a pure vacuum or filled with conducting plasma, and 
whether or not we include viscosity. In the first case (vacuum 
exterior, no viscosity) our treatment is entirely analytic, while 
in the other two the analyses are numerical. In Section 5 we 
restore dimensions and evaluate our key results for the case of a 
differentially rotating remnant of a binary neutron star merger. 
Here we summarize the implications of our calculations to fu- 
ture observations and numerical simulations. Appendicies A 
and C summarize details of our numerical calculation, while 
Appendix B provides a proof that the lowest energy state for 
our adopted rotating stellar models is the state of uniform rota- 
tion, for fixed angular momentum. 

2. BASIC MODEL 

As discussed in BSS, differential rotation in a spinning star 
twists up a frozen-in, poloidal magnetic field, creating a very 
strong toroidal field. This process will generate Alfven waves, 
which can redistribute and even carry off angular momentum 
from a star on timescales less than ~ 100s for poloidal fields 
greater than ~ 10 I2 G. Shear viscosity also redistributes angular 
momentum. However, molecular viscosity in neutron star mat- 
ter operates on a typical timescale of ~ 10 9 s in a star of ~ 10 9 K, 
so it alone is much less effective in bringing the star into uni- 
form rotation than magnetic braking, unless the initial magnetic 
field is particularly weak. Turbulent viscosity, if it arises, can 
act more quickly. 

To track the evolution of differential rotation and follow the 
competition between magnetic braking and viscous damping, 
we model a spinning star as an infinite, axisymmetric cylinder. 
We take the gas interior to the star to be homogeneous in den- 
sity and incompressible, which is a crude but not unreasonable 
approximation for a neutron star governed by a stiff nuclear 
equation state. We adopt cylindrical coordinates (r,<j>,z), with 
the z— axis aligned with the rotation axis of the star. We take 
the magnetic field to have components only in the r and <\> direc- 
tions. We consider two cases: one in which the exterior region 
of the star is pure vacuum and another in which it is a diffuse 
plasma at a constant density. We treat the fluid everywhere to 
be perfectly conducting and allow for the presence of shear vis- 
cosity, which we take to be constant. 

A similar geometry was adopted previously to study mag- 



netic braking of a rotating, cylindrical gas cloud immersed in 
the interstellar medium (Mouschouvias & Paleologou 1979). 
However, in that analysis it proved sufficent to evolve the mag- 
netic field only in the exterior, treating the cloud as a rigid rota- 
tor and ignoring viscosity altogether. Here we are specifically 
interested in the evolution of the field inside the star and learn- 
ing exactly how the field and viscosity combine to alter the 
interior angular momentum profile of a differentially rotating 
configuration. 

3. BASIC EQUATIONS 

The fundamental equations for our incompressible, perfectly 
conducting MHD fluid are the equation of continuity 



V-v = 0, 



(1) 



where v is the velocity, and the magnetic Navier-Stokes equa- 
tion 



dv 



+ (v V)v = 

1 , (V x B) x B 

VP - V<f> + jA7 2 v + , 

p 4np 



(2) 



where p is the density, v is the viscosity, P is the pressure, B is 
the magnetic field and $ is the Newtonian gravitational poten- 
tial, which satisfies Poisson's equation, 

V 2 $ = 4np. (3) 

The magnetic field B satisfies Maxwell's constraint equation 

V-B = 0, (4) 

as well as the flux-freezing equation 

<9B 

— = V x (v x B). (5) 

Solving equation ([!]) for our cylindrical star gives the radial 
component of the velocity as v r = Rvo(t)/r where R is the stel- 
lar radius. Demanding the flow be regular at the origin requires 
that the function Vo(t) vanish identically, implying that only ro- 
tational motion can take place, with an azimuthal velocity given 
by = r£l(t,r). Here Q, is the angular velocity. Solving equa- 
tion (Q) for the magnetic field together with the flux-freezing 
condition (Q) requires that the radial component of the magnetic 
field be independent of time and given by 



B r = 



B a R 



t >0, 



(6) 



where Bq is the value of the field at the stellar surface at t = 0. 
Although the magnetic field given by equation (||) exhibits a 
static line singularity along the axis at r = 0, it does not drive 
singular behavior in the fluid velocity or nonradial magnetic 
field. As these quantities, which are the main focus here, re- 
main finite and evolve in a physically reasonable fashion, the 
line singularity poses no difficulty and requires no special treat- 
ment. 

Evaluating the azimuthal components of equations (Q) and 
(||), the system then reduces to two coupled partial differential 
equations for £1 and Z?^, 



Oil _ RB d 
dt 4irpr 3 dr 



rB ri 



v d 
r 3 dr 



(7) 
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and 



dt 



on 

or 



(8) 



Equations (0) and (|^) are the main dynamical equations to be 
solved. Once and B^ are determined, the radial component 
of the Navier- Stokes equation, 



dP 



9$ 



idBl 



dr ^ r ^ dr 4-tt \ 2 dr r 



(9) 



can be solved for the pressure as a function of r and t by a 
simple integration over r. Though straightforward, obtaining 
the resulting pressure, which is the profile necessary to enforce 
hydrostatic equilibrium in the star in the radial direction at all 
times, is not required to obtain the evolution of the rotation pro- 
file or magnetic field and is therefore of little interest here. Sim- 
ilarly, the equation of heat transfer is not needed in solving the 
problem of incompressible flow, as we are not interested here 
in the temperature distribution. 

3.1. Initial Conditions and Boundary Values 

We assume that no azimuthal component of the field B$ is 
present initially, 

B*(0,r) = 0, r>0, (10) 

recalling that B, is always given by equation (||). We are thus 
interested in the situation where the azimuthal field is created 
entirely by the differential rotation of the fluid, which bends the 
frozen-in, initially pure radial field lines in the azimuthal direc- 
tion. 

We solve equation ^ for fl subject to the condition of regu- 
larity at the origin, 



dfl(t,0) 
dr 



= 0, t >0, 



which, by equations (|^) and (|lC|), implies 
dB4t,Q) 



dt 



= = B^(t,0), t>0. 



(11) 



(12) 



The boundary conditions at the stellar surface depend on the 
physical situation being investigated. In the case of an exterior 
vacuum, no azimuthal magnetic field can be carried into the 
region outside the star. This fact, together with equation ([$]), 
implies 

dfl(t,R) 

BAtjR) = = — , t >0 (vacuum exterior). (13) 

or 

In the case where the exterior region of the star contains a 
conducting homogeneous fluid at density p ex , part of the az- 
imuthal field is transmitted across the stellar surface in the form 
of an outgoing Alfven wave, and part is reflected via an ingo- 
ing Alfven wave. We construct an appropriate wave boundary 
condition accounting for partial transmission and reflection at 
the interface between the interior star and the exterior plasma 
according to 

dB^{t,R) dB^t,R) (1+ft) 

H ^ va — — — = (plasma exterior), (14) 



dt 



dr 



where 1Z is the reflection coefficient for the azimuthal field am 
plitude, 

0Wp) 1/2 -i 



0Wp) 1/2 +i' 



(15) 



and va = Bo/(47rp) 1 / 2 is the Alfven speed just inside the stellar 
surface (see Appendix A). The corresponding boundary condi- 
tion for 51 is 

dn(t,R) dsi(t,R) i q+g) 

H — — — = (plasma exterior). (16) 



dr 



dt v A (l-Tl) 



Note that when p ex = 0, we have 1Z = -1, so equations ( |14f) and 
( |i"oj ) recover the vacuum exterior boundary conditions (jT3|), re- 
calling equation (|To|). By contrast, when p ex = p, and the ex- 
terior plasma is a smooth continuation of the interior star, we 
have TZ = 0, and equations (n~4b and (16) impose standard out- 



going (Alfven) wave boundary conditions on these transverse 
field and velocity components. 

The entire evolution of the system is driven by differential 
rotation. We take the initial rotation profile of the star to have 
the following nonuniform, momentarily stationary form: 

fi(0,r) = ^ [l + cos(7rr 2 /fl 2 )], < r < R, (17) 



where 



aO(0,r) 
dt 



0, 0<r<R. 



(18) 



We note that our adopted initial angular velocity profile is con- 
sistent with the boundary conditions at the center and surface 
of the star (see equations ( fi"T| ) and (fl3|)), as required. 

3.2. Conserved Energy and Angular Momentum 

The MHD equations (Q) and admit two nontrivial inte- 
grals of the motion, one expressing conservation of energy and 
the other conservation of angular momentum of the star. Energy 
conservation is given by the integral 



£ rot (0) = ErJt) + £ mag (f) + £ vis (?) + £ Poyn (f), 



(19) 



where E m (t) is the rotational kinetic energy of the matter, 
£mag(f ) is the azimuthal magnetic energy, £ v is(0 is the internal 
(thermal) energy generated by viscous dissipation and Ep oyn (t ) 
is the energy carried off at the stellar surface by the Poynting 
vector S = cE x B/47T, where E = -v x B/c: 

E rot (t) = J d\( P n 2 (t,r)r 2 /2),E mag (t) = J tfxXtffts)/**) 



E Yis (t)= / dt£ vis (t),E Poya (t)= / dt£ Poyn (t). (20) 
Jo Jo 

The rates of viscous dissipation £ v ^(t) and Poynting energy loss 
<£p yn(f ) appearing in equation (|o|) are given by 



t\ .jn = I d 3 x 

Bf^By^lR 



dfl 



(21) 



£poyn(0 



4n 



A. 



The volume elements in equations ([Z0|) and (g2|) are evalu- 
ated as d^x = 2irLrdr, with the volume integration performed 
throughout the interior of the cylindrical star over a (arbitrary) 
height L in the z-direction. The Poynting flux is measured at 



4 



the surface of this region, whose area is A = 2ttLR. Dividing all 
quantities by L then gives energies per unit length. 
Conservation of angular momentum is expressed as 2 



/rot(O) = /rot(0 + /mag(0, 



(22) 



where J mt is the rotational angular momentum of the matter and 
7 mag is the integrated torque exerted by the Maxwell stress at the 
surface, 

/rotW = J d 3 x(pQ(t,r)r 2 ), 7 mag (f) = J dtM ". (23) 
The torque TV is given by 



Af = 



47T 



A 



(24) 



We note that, consistent with the nonrelativistic MHD approxi- 
mation, the electric field energy E 2 /8tt is not included in equa- 
tion (19) and the angular momentum of the electromagnetic 
field S<f,/c 2 is not included in equation (22) (Landau, Lifshitz 
& Pitaevskii 1984). 

The motivation for the monitoring the conservation equations 
during the evolution is twofold: physically, evaluating the indi- 
vidual terms enables us to track how the initial rotational energy 
and angular momentum in the fluid are transformed, dissipated 
and/or transported away; computationally, monitoring how well 
the conservation equations are satisfied provides a check on the 
numerical integration scheme. 

3.3. Nondimensional Formulation 

We introduce nondimensional variables to facilitate the inte- 
gration and to help identify the scaling behavior of our solution 
with respect to our arbitrary choices of input parameters. We 
define nondimensional quantities according to 

r=r/R, i=2t/(R/v A ), u = 4u/(v A R) (25) 

Q = n/Q , B0 = B0/((47rp) 1 / 2 (n o i?)) 

E = E/ ((irR 2 Lp)(n 2 R 2 /2)) , J = J/ ((nR 2 Lp)(n () R 2 )) . 

We work with nondimensional variables in all subsequent equa- 
tions, but for simplicity we omit placing carets O on the vari- 
ables. In nondimensional variables, the coupled evolution equa- 
tions ([7]) and (||) become 



an 1 d 



8 



dt 



dr 4 V dr 2 



and 



OBcj, _ an 

~W r dr 2 ' 



(26) 



(27) 



where the nondimensional energy variables now appearing in 
that equation are 



E mt (t) = fudr 2 n 2 (t,r)r 2 , £ mag (f) = ^dr 2 B 2 
E^{t) = v&dt [j> 2 (fl)V 
£poyn(0 = _ 2 JydtiB^Q),-^ . 

Multiplying equation d27| ) by r 2 dr 2 and integrating over the star 
yields 2 angular momentum conservation equation (p2[), where 
the nondimensional angular momentum terms are given by 

/«*(«)= / dr 2 fl(t,r)r 2 , 7 mag (f) = - f dt{B^) r=l . (29) 
Jo Jo 

The coupled system (|6|) and (|7|) is purely hyperbolic in 
the absence of viscosity and leads to a simple wave equation. 
The inclusion of viscosity gives the system a mixed hyperbolic- 
parabolic form, since viscosity contributes diffusive damping to 
the evolution. 

4. SOLUTIONS 

We solve the evolution equations for three different cases. In 
case A, we treat a viscous-free star with a vacuum exterior. In 
case B, we consider the effect of viscosity. In case C, we again 
ignore viscosity but assume that the exterior region consists of 
a perfectly conducting plasma of lower density than the central 
star. By examing these three cases separately we are able to dis- 
tinguish the different physical processes by which the rotation 
profile in the star is altered. 

4. 1 . Case A: v = with vacuum exterior 

Here we treat the simplest case in which the fluid is assumed 
to have no viscosity and the region exterior to the star is a pure 
vacuum. Setting c = Owe take the partial derivative of equation 
( [26| ) with respect to time and use equation (27 ) to eliminate B$, 
thereby obtaining a simple wave equation for fi, 

d2Q ^r,on\ (30) 



(28) 



dt 2 



i__d_ 

r 2 dr 2 



dr 



This equation must be solved subject to the boundary condi- 
tions (|nj) and @ 



dSl(t,Q) 



= = 



dCl(t,l) 



dr dr 
and initial conditions (17) and ( |I"8| ) 



(1(0,0 = ifi+cosc^)], m0 > r) 



dt 







(31) 



(32) 



We recognize equation (30) as Oil = expressed in cylindri- 
cal coordinates, with cylindrical radial coordinate £ = r 2 and 
no dependence on <fi or z. In these coordinates, the equation 
is identical to the one that would arise when determining the 
amplitude of the oscillations of a circular, axisymmetric drum- 
head, undamped at the outer edge (see, e.g, Mathews & Walker 
1970 for an analysis of a circular drumhead that is clamped at 
the outer edge). We solve the equation analytically by separa- 
tion of variables, obtaining an expansion over Bessel functions 
according to 



When expressed in nondimensional units, the initial value equa- 
tions and boundary conditions take on the same appearance as 
in equations © - (§J|), provided we set v'a = 1 = R wher- 
ever these quantities appear. Multiplying equation (^) by 
flr 2 dr 2 , integrating over the radius of the star and making use 
of use equation (|7j) yields energy conservation equation (|T9|), 

2 In deriving both equations ( |l9| ) and ( |2^ ) we have set the viscous surface terms proportional to vdVljdr equal to zero, which is appropriate for all of the cases 
treated below. 



fi(f,r) = ^2 = A n Ja(k n r 2 )cos(k n t), 



(33) 
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where 



A n 



2^dr 2 r 2 J (k n r 2 )n(0,r) 



(34) 



and where Q(0,r) is set by equation ([32|). In the above expan- 
sion, the quantities {k n } , n = 1 , 2 , . . . are the zeroes of the deriva- 
tive of Jo(x), i.e., they satisfy dJo(k n )/dx = -Ji(k n ) = and are 
tabulated (Abramowitz and Stegun (1972), Table 9.5, p. 409, 
remembering to add k\ = 0). The azimuthal field is obtained by 
substituting equation ( ^3| in to equation (p7|), integrating over 
time and using equation (|lO|) to set the initial value. The result 
is 



B <t> = -r^A n J\(k n r 2 )sm(k„t). 

n=l 



(35) 



We evaluate the coefficients A„ up to n = 20 by quadrature 
and the resulting standing wave solutions for the angular ve- 
locity and azimuthal magnetic field are shown in Figs. 1-3. 
In Fig. 1 we track the evolution of the various energies with 
time and observe how the total conserved energy oscillates be- 
tween rotational kinetic energy of the fluid and energy of the az- 
imuthal magnetic field. The oscillations proceed forever with- 
out any dissipation or energy loss. Angular momentum, which 
is carried entirely by the fluid, is also strictly conserved. We 
plot the angular velocity profile at critical phases during an os- 
cillation period {P = 1 .64) in Fig. 2 and we plot the correspond- 
ing azimuthal field (a standing Alfven wave) in Fig. 3. 

The series solutions ( f33| ) and ( J35] ) converge rapidly: the first 
few values of A n are A x = 0.2974, A 2 = 0.7192, A 3 = -0.0198 
and A\ = 0.0044 . We already obtain a reasonable approxima- 
tion by using the first two coefficients alone, together with the 
quantities k\ = and k 2 = 3.8317, to get 

ft(f,r)«0.2974 + (0.7192)/ (3.8317r 2 )cos(3.8317f), (36) 



and 



-(0.7192)^(3.8317 r 2 )sin(3.8317 t). 



(37) 
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FIG. 1. — Energy evolution for a differentially rotating star with zero vis- 
cosity and a vacuum exterior. The dotted line shows £Vot and the dashed line 
shows £ m ag- All energies are normalized to E rot (0). The sum of the energies 
is conserved and remains equal to the initial rotational energy E mt (0), which is 
plotted as the solid horizontal line at Energy = 1 . Time is in nondimensional 
units. 
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FIG. 2. — The angular velocity profile at critical phases during the first oscil- 
lation cycle for the differentially rotating star described in Fig. 1. Curves are 
labeled by the value of time; the standing wave pattern oscillates with a period 
P = 1.64. All quantities are plotted in nondimensional units. 
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FIG. 3. — The azimuthal magnetic field profile at critical phases during the 
first oscillation cycle for the differentially rotating star described in Fig. 1. 
Curves are labeled by the value of time; the standing Alfven wave pattern 
repeats with a period P = 1 .64. All quantities are plotted in nondimensional 
units. 



From the above approximation we derive the oscillation period 
to be P m 2n/k2 = 1.64, consistent with the behavior shown in 
Figs. 1-3. From equation ( ]36| ) we calculate f2 U nif, the value 
of the angular velocity at the phases of uniform rotation, to 
be O un if = A\ = 0.2974, consistent with the horizontal lines in 
Fig. 2. This value of 57 can also be found by invoking angular 
momentum conservation: setting 7 ro t(0), the total (initial) angu- 
lar momentum of our differentially rotating star, obtained using 
equations (^9|) and (f32|), to the angular momentum of a rigidly 
rotating configuration with the same density profile, we find 



n 



unif " 



2tt 2 



: 0.2974 



(38) 



At these phases achieves its maximum amplitude. The en- 
ergy at the uniformly rotating phases is distributed more evenly 
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between rotational kinetic and azimuthal magnetic field en- 
ergy. In our model, the rotational energy falls to a fraction 
£rot/£rot(0) = 2(7r 2 -4) 2 /(7r 2 (37r 2 - 16)) = 0.513, while the mag- 
netic energy increases to E m!lg /E ml (Q) = 1 -E lol /E mt (0) = 0.487 
at these phases. In the absence of dissipation, the star oscillates 
about the uniformly rotating state indefinitely. 

The presence of internal energy dissipation (e.g., viscosity), 
however small, will drive a differentially rotating star to a per- 
manent state of uniform rotation, since uniform rotation is the 
lowest energy state at fixed angular momentum (see Appendix 
B for a proof). In this situation, the dissipated energy ultimately 
goes into heat, and the azimuthal magnetic field decays away. 
We will study this behavior in Case B, below. However, even 
in the absence of dissipation, the presence of a toroidal (here, 
radial) magnetic field will brake the differential rotation on the 
time it would take an Alfven wave, traveling at the local Alfven 
speed VaM = vp,(R)R/r, to cross the radius of the star. If we 
define ?a = R/va, then in the model analyzed here the Alfven 
wave crossing time is exactly = 1 in nondimensional units. 

The scaling behavior of our solution, revealed by our nondi- 
mensional formalism, exhibits several significant features. We 
find that the values of the amplitudes of all evolved quantities 
shown in Figs. 1 -3 are entirely independent of the magnitude 
of the initial radial seed field Bo given by equation (^). Only 
the timescale of the evolution depends on the strength of the 
seed field, and it is proportional to the Alfven time associated 
with that field. Specifically, no matter what the strength of the 
initial radial field, the azimuthal magnetic field will grow to the 
same high value sufficient to brake the differential motion and 
drive the star to oscillate about the state of uniform rotation. In 
this state, the energy in the azimuthal field will always grow 
to make up the difference between the rotational energy in the 
initial, differentially rotating configuration and the uniformly 
rotating configuration, independent of the strength of the seed 
field. By the same token, the evolution timescale is independent 
of the degree of differential rotation or the rotation period. 
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FIG. 4. — Energy evolution for a differentially rotating star with internal vis- 
cosity v — 0.2 and a vacuum exterior. The dotted line shows E Iot , the dashed 
line shows E mag , and the long-dash-dotted line shows E v [ s . All energies are 
normalized to E ro t(0). The sum of the energies, plotted as a short-dash-dotted 
line, is conserved and remains equal to the initial rotational energy £ rol (0), 
plotted as the solid horizontal line at Energy = 1 . Time is in nondimensional 
units. 



4.2. Case B: Nonzero viscosity with vacuum exterior 

We now consider the effect of viscous dissipation by setting 
v = 0.2 in equation (p6|). This intermediate value is chosen in 
order that viscosity in the star be sufficiently large for dissi- 
pation to become evident in a few Alfven timescales, but suf- 
ficiently small so that viscous damping of differential rotation 
does not completely suppress the growth of the azimuthal mag- 
netic field. We solve the coupled evolution equations ( |26| ) and 
(p7[), subject to the same vacuum exterior boundary and ini- 
tial conditions used in Case A: equations ( |lO| ) - (13) for B^, 
together with (|Tj) and (|2])) for f2. We solve these equations 
numerically by finite-differencing (see Appendix C). To better 
treat the multiple timescales (e.g., Alfven vs. viscous) and to 
follow the long-term damping, we employ an implicit scheme 
and thereby eliminate a Courant stability bound on our numeri- 
cal timestep. We test our numerical integrations by reproducing 
the analytic solution derived in Case A for v = and by check- 
ing that energy and angular momentum are conserved accord- 
ing to equations (|19|) and @. 

The results of our numerical integrations are summarized 
in Figs. 4-6. In Fig. 4 we track the evolution of the 
various energies with time. We observe how the oscilla- 
tions of the rotational kinetic and azimuthal magnetic field 
energies are now damped by viscosity. From the evolu- 
tion equation (|26|), we may define the characteristic viscous 
dissipation timescale according to t v = v~ l in nondimen- 
sional units, which is equivalent to t v = R 2 /(8z/) in physi- 
cal units (see equation (E5b). For the example considered 
here, the ratio of viscous to Alfven time is t u /tA = l/(2f) = 
2.5. The numerical results are consistent with this ratio. 
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FIG. 5. — The angular velocity profile at critical phases during the first three 
oscillation cycles for the differentially rotating star described in Fig. 4. Curves 
are labeled by the value of time. The wave pattern undergoes damped oscil- 
lations with a period P = 1.64. All quantities are plotted in nondimensional 
units. 

As time progresses, the rotation becomes uniform (see Fig. 5), 

the azimuthal field decays back to zero (see Fig. 6), and the 
rotational energy difference bet ween the initial differentially 
rotating and the final uniformly rotating star is dissipated as 
internal heat by the viscosity. Angular momentum is again con- 
served in the star, as verified by the numerical results. Fig. 4 
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shows that at late times the final rotational energy approaches 
E rot /E mt (Q) = 0.513, the value already calculated in Case A for 
uniform rotation for the given (conserved) angular momentum 
of our star. At the same time, the en ergy in viscous dissipation 
increases to E v i s /E rot (Q) = 1 -E rot /E mt (0) = 0.487. Similarly , 
Fig. 5 shows that undergoes damped oscillations about the 
asymptotic, uniform value f2 un if give n by equation (|3"s|). 

The asymptotic stationary state of the star, which consists of 
uniform rotation, zero azimuthal magnetic field, and apprecia- 
ble internal heat content is uniquely determined by conserva- 
tion of total energy and angular momentum: it does not depend 
on the strength of the assumed viscosity, only that some vis- 
cosity is present, however small. The magnitude of the vis- 
cosity does determine the damping timescale for differential 
rotation and Alfven wave oscillations. For a sufficiently high 
viscosity v ^> 1, differential rotation is damped without driving 
Alfven waves throughout the star. However, the final state is 
unchanged. 
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lar momentum to the exterior plasma. Fig. 9 shows that the an- 
gular velocity profile becomes increasingly uniform with time 
and that the star ends up spinning in the opposite sense from its 
rotation at t = 0. By the time the star settles into uniform rota- 
tion, the interior azimuthal field dies away, as shown in Fig. 10. 
Figs. 7 and 8 confirm that total energy and angular momentum 
are conserved, in compliance with equations ( fl^ ) and I 
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FIG. 7. — Energy evolution for a differentially rotating, zero-viscosity star 
embedded in an ambient plasma where pex/p = 0.2. The dotted line shows 
E I0 u the short-dashed line shows E ma g, and the long-dashed line shows the 
outgoing integrated Poynting energy flux £p oyn . All energies are normalized 
to £V ot (0). The sum of the energies, indicated by the dot-dashed line, is con- 
served and remains equal to the initial rotational energy £ ro t(0), plotted as the 
solid horizontal line at Energy = 1. Time is in nondimensional units. 
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FIG. 6. — The azimuthal magnetic field profile at critical phases during 
the first three oscillation cycles for the differentially rotating star described in 
Fig. 4. Curves are labeled by the value of time. The Alfven wave pattern un- 
dergoes damped oscillations with a period P = 1.64. All quantities are plotted 
in nondimensional units. 



4.3. Case C: v = with plasma exterior 

Here we treat a differentially rotating star that has no viscos- 
ity but is surrounded by an infinite, homogeneous medium of 
perfectly conductingplasma. We solve the coupled evolution 
equations ( ^6| ) and Q27| ) for the same boundary and initial con- 
ditions as in Case A, except at the surface, where we now use 
equation (JIJ) in place of ( fl3| ) to account for the partial trans- 
mission and reflection of the Alfven wave at the star - plasma 
interface. We set p ex /p = 0.2. As in Case B, we again solve the 
equations by finite differencing (see Appendix C). 

Our results are summarized in Figs. 7 - 10, where we see 
how the loss of energy and angular momentum to the exterior 
via the transmission of toroidal Alfven waves drives the star to 
a slower, rigidly rotating state. In Fig. 7 we observe how the 
oscillatory exchange of energy between rotation and toroidal 
magnetic field is damped by the outgoing Poynting energy flux 
at the surface. Fig. 8 shows how the torque exerted on the star 
by the magnetic field at the surface causes the star to lose angu- 
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FIG. 8. — Angular momentum evolution for the differentially rotating star 
depicted in Fig. 7. The dotted curve shows the angular momentum of the star, 
while the dashed curve shows the integrated torque exerted by the Maxwell 
stress at the stellar surface. All contributions are normalized to /rot(0). The 
sum of all angular momenta, a dot-dashed line, is conserved and remains equal 
to 7 ro i, plotted as the solid horizontal line at Ang Mom = 1. 
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The existence of an ambient plasma is sufficient to drive the 
spin of the star to low values; internal dissipation is not neces- 
sary. In this model, the magnitude of the final (retrograde) spin 
of the star, ss -0.1, is insensitive to the ratio of the densities 
of the homogeneous ambient plasma to the homogeneous stel- 
lar interior, p ex / p. This ratio influences the timescale of rotation 
damping only. A crude estimate suggests that the net damping 
timescale is roughly ?A,ex ~ £a(/V A*) 1 » f° r small values of 
p ex / p. This estimate accounts for the fact that only a fraction 
« 4(p ex /p) 1//2 of the outgoing Poynting energy flux carried by 
the Alfven wave actually escapes from the stellar surface into 
the ambient gas during each oscillation cycle, while the rest is 
reflected back inward (see equation [ A4]). 
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FIG. 9. — The angular velocity profile at critical phases during the first three 
oscillation cycles for the differentially rotating star described in Fig. 7. Curves 
are labeled by the value of time. The wave pattern undergoes damped oscil- 
lations with a period P = 1.64. The star settles into retrograde motion with 
angular velocity f2 R2 — 0. 1 , indicated by the dashed line. All quantities are 
plotted in nondimensional units. 



5. DISCUSSION 

The presence of a seed poloidal magnetic field in a differen- 
tially rotating star may be sufficient to drive the star to a state 
of uniform rotation. Differential rotation generates toroidal 
Alfven waves through the star and these waves can convert an 
appreciable fraction of the kinetic energy in differential mo- 
tion into magnetic field energy. Even in the absence of inter- 
nal dissipation or Poynting energy loss to an external plasma, 
the conversion of rotational to magnetic energy by these waves 
can remove a portion of the star's rotational support against 
gravity, which can upset the equilibrium balance in a station- 
ary configuration. In the case of a hypermassive neutron star 
supported against collapse by differential rotation, the gener- 
ation of Alfven waves can lead to catastrophic collapse to a 
black hole, possibly accompanied by mass loss. Such an out- 
come may characterize the final fate of neutron stars formed in 
binary mergers (BSS). 

Our model calculations help to identify at least five dis- 
tinct timescales that govern the evolution of a a differentially 
rotating, equilibrium neutron star containing a seed poloidal 
magnetic field. We evaluate them here for neutron star pa- 
rameters consistent with the stable remnant found in the rel- 
ativisitic simulations of binary mergers by Shibata and Uryu 
(2000). These parameters are comparable to the values of the 
hypermassive, differentially rotating equilibrium model con- 
structed and evolved by BSS, who found that their configura- 
tion was also dynamically stable. All of these models employ 
a poly tropic equation of state with poly tropic index n = 1 , for 
which there is scale freedom to choose the physical mass. If 
we take the mass to be M w 2 x 1.5M0 3M0, a reason- 
able value for a binary remnant, then the equatorial radius of 
the star is R « 5M w 20km and the central angular velocity is 
flo ~ 0.3/M w2x 10 4 Hz, a factor of three times higher than 
the angular velocity at the equator. 

The dynamical timescale associated with the star is given by 
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FIG. 10. — The azimuthal magnetic field profile at critical phases during 
the first three oscillation cycles for the differentially rotating star described in 
Fig. 7. Curves are labeled by the value of time. The Alfven wave pattern un- 
dergoes damped oscillations with a period P = 1.64. All quantities are plotted 
in nondimensional units. 



the time it takes the binary remnant to achieve equilibrium fol- 
lowing coalescence. This is also the time that it takes the star, 
once it is driven out-of-equilibrium by magnetic braking, to 
undergo collapse. The central rotation period of the remnant is 
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while the period at the equator is about three times longer. 
The timescale for magnetic braking of differential rotation by 
Alfven waves is given by 



R 

VA 



10- 
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(41) 



On this timescale the angular velocity profile in the star is sig- 
nifcantly altered, and this change can drive the star far out of 
equilibrium whenever the initial rotation supplies an apprecia- 
ble fraction of the support against gravity. In the case of a hy- 
permassive neutron star, which depends on rapid differential 
rotation for hydrostatic support, magnetic braking can lead to 
catastrophic collapse to a Kerr black hole, possibly accompa- 
nied by some mass and field ejection. A hypermassive star is 
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a likely outcome of binary neutron star coalescence. The re- 
sulting delayed collapse following merger may generate a brief 
secondary burst of gravitational waves lasting a time fdyn- If the 
time of the coalescence can be inferred from the initial wave 
burst signal, then the measurement of this delay should provide 
an estimate of the strength of the poloidal magnetic field in the 
interior of the merged neutron star (BSS). 

Neutron stars below the supramassive limit do not require 
differential rotation for hydrostatic support. Such stars have 
masses at most 20% larger than the nonrotating, spherical TOV 
limit (see, e.g, Butterworth & Ipser (1975); Friedman, Ipser & 
Parker (1986); Komatsu, Eriguchi and Hachisu (1989); Cook, 
Shapiro & Teukolsky 1992,1994; Salgado et al. 1994). 3 If 
such a star is formed with differential rotation, magnetic brak- 
ing will still drive Alfven waves in the star and alter the rotation 
profile on the Alfven time. While the star may undergo dynami- 
cal oscillations, it will not collapse. Instead, viscous dissipation 
ultimately will drive the star to a new, uniformly rotating equi- 
librium state on a timescale 
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(42) 
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where v = 347 ( o 9 / 4 r- 2 cmV 1 (Cutler & Lindblom 1987). 

If a differentially rotating neutron star is immersed in an am- 
bient plasma at formation, Alfven waves generated in the inte- 
rior will propagate into the external medium, carryin g of f en- 
ergy and angular momentum. According to equation ( A4), ev- 



erytime a wave propagates to the stellar surface it will transmit 
a fraction w 4(p ex /p) 1//2 , (p ex / p) <C 1, of its energy to the sur- 
rounding plasma, causing the star to spindown in an e-folding 
time 
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(cf. Spitzer 1978, eqn 13.54 and discussion therein for an alter- 
native argument leading to an analogous result). Determining 
the final spin rate and fate of a newly formed neutron star may 
thus depend on the nature of its immediate surroundings. 

Our naive model calculations do not account for convective 
instabilities (Pons et al. 1999), which can drive the seed mag- 
netic field to high values greatly exceeding 10 12 G (Duncan & 
Thompson 1992), or other possible MHD instabilities, which 



may also contribute to the redistribution of angular momentum 
(Balbus & Hawley 1994; Spruit 1999b). Our main goal was 
to show that even in the simplest laminar description, a differ- 
entially rotating neutron star is a transient phenomenon, since 
magnetic braking and viscosity inevitably will bring a differ- 
entially rotating star into uniform rotation (although not be- 
fore many dynamical timescales have elapsed, according to our 
timescale estimates above). Thus all radio pulsars are likely to 
be uniformly rotating. However, differential rotation may char- 
acterize a nascent neutron star formed in a supernova, following 
fallback, or in a merged binary. The braking of this motion may 
have important consequences for gravitational wave signals and 
gamma-ray bursts. 

More realistic evolutionary calculations of magnetic brak- 
ing in neutron stars should clarify some of the above issues. 
One computational subtlety arises from the inequality fdyn *C tA, 
which holds for seed magnetic fields of small or moderate 
strength (see equations [ f39[ ] and [^TJ). In this regime it may 
prove too taxing to a relativistic MHD code to evolve a differ- 
entially rotating star for the required length of time for magnetic 
braking to take effect. An obvious solution would be to artifi- 
cially amplify the initial magnetic field strength, increasing the 
ratio fdyn/fA to a computationally manageable size (while still 
keeping it much less than unity), and then scale the results for 
smaller ratios. However, if the energy in the seed magnetic field 
is an appreciable fraction of the gravitational potential energy, 
care must be taken to incorporate the magnetic field stresses 
when constructing the initial hydrostatic equilibrium configu- 
ration. Another approach might be to use implicit differencing 
in the calculation to avoid the Courant stability constraint on 
the evolution timestep. In addition, one might be able treat part 
of the evolution in the quasistatic approximation, as in a typical 
stellar evolution code, up to the moment that stable equilibrium 
can no longer be sustained. Specifically, one might evolve the 
magnetic field and the angular velocity on the Alfven timescale 
via the MHD equations, but assume hydrostatic equilibrium is 
maintained on the dynamical timescale in order to determine 
the density and pressure profiles at each instant. Our idealized 
incompressible cylindrical star model, which remains radially 
static, trivially conforms to this behavior and thereby avoids the 
mismatch between the dynamical and Alfven timescales. (An 
additional motivation for treating this idealized problem was to 
highlight this feature of the analysis.) We hope to tackle the 
more general problem in future computational investigations. 

It is a pleasure to thank T. Baumgarte, C. Gammie and M. 
Shibata for helpful discussions. This work was supported by 
NSF Grants AST 96-18524 and NASA Grants NAG 5-7152 and 
NAG 5-8418 at the University of Illinois at Urbana-Champaign. 



APPENDIX 

PARTIALLY TRANSMITTED/REFLECTED WAVE OUTER BOUNDARY CONDITION 

If the surface of a star is the interface between a highly conducting interior plasma of density p and a highly conducting exterior 
plasma of density p ex , an Alfven wave generated in the interior progating to the surface suffers partial reflection and transmission 
(Roberts 1967). The reflection and transmission coefficients for the transverse magnetic field amplitude are 

* - 'M^r (Al) 

(Pex/p) 1/2 +l 

3 The most recent equations of state place the maximum gravitational mass of a spherical star in the range 1.8-2.2Mq (Akmal, Pandharipande & Ravenhall 1998), 
but observations of SN 1987A suggest that it could even be as low as 1.5Mq (Brown & Bethe 1994) 
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and 



T = 2{P :'^ 2 (A2) 



= ( S PgX / P !l/2 , ! ) (A3, 



The relection and transmission coefficients for the Poynting energy flux are 

>ex/p) 1/2 +l 

and 

Tp „„, . - ^f. . (A4) 

((fc,/p)" ! +l) 

Note that 7\Lp oyn + 7p oyn = 1, in compliance with energy conservation. 

We seek to impose an approximate boundary condition on the transverse magnetic field just inside the surface. Decomposing the 
wave amplitude of the field into its outgoing incident and ingoing reflected components, we may write 

B^t , r) = f(t - r/v A ) + Kf(t + r/v A ) (A5) 



where va is the Alfven propagation speed at the surface. In writing equation (A.5) we approximate the wave to be planar, which is 
strictly true only for wavelengths Aa <C R- Differentiating equation (A5) yields 

1 gg^ v A 

J dt (1 + ft) dr (I-IZ) K ' 

where the prime denotes a total derivative. We may then recast equation ( |A6| ) as 

dB,p(t,R) i 8B4t,R) (1 + K) 

+ 5 Va 7i — t?\ =U ' ( A/ - ) 

at or 

which provides the desired boundary condition. 

UNIFORM ROTATION AS THE LOWEST ENERGY CONFIGURATION 

Here we prove that for fixed angular momentum, the angular velocity profile that gives a stationary star of lowest energy is one 
of uniform rotation. To find the lowest energy state, set the azimuthal magnetic field energy to zero and take the total energy of the 
configuration to be rotational kinetic energy, 

E= [ dr 2 E(n,r 2 )= [ dr 2 fl 2 r 2 , (Bl) 
Jo Jo 

The conserved angular momentum is given by 

J= [ dr 2 J(fl,r 2 )= [ dr 2 flr 2 . (B2) 



The lowest energy configuration is found by varying the energy subject to the constraint of fixed angular momentum, i.e., by varying 
the integral 

1 = E + XJ (B3) 
where A is a (constant) Lagrange multiplier. Setting the variation equal to zero yields the Euler-Lagrange equation for the functional 

C(n, fl',r 2 ) = E(Q,r 2 ) + XJ(Q, r 2 ). (B4) 

The Euler-Lagrange equation gives 

— = (2ft+A)r 2 = 0, (B5) 

or 

ft = —— = constant, (B6) 

which is the desired result. 
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FINITE DIFFERENCE EQUATIONS 

To finite difference the coupled evolution equations (^) and ( |27| ) we introduce a uniform radial grid r,-, i= 1,2, ... ,i max , which 
extends from r\ = to r, nrat = 1 . The angular velocity f2 is defined at the midpoints of the radial zones while the magnetic field B is 
defined on the zone boundaries. We adopt an implicit evolution scheme in which the differencing is first order in time but second 
order in space. Care is taken when finite differencing the radial gradients so that they reflect the correct regularity behavior at the 
origin, where B ~ a\r, while £1 ~ a^ + a^r 1 ; here ao, a\ and «2 are constants. 

To evolve f2 for a timestep At = t n+ \ - t n we difference equation (Eq) according to 



QM+l _ Qll 

"i+1/2 "f+1/2 

At 



1 


(n +l B';:l- n B1 +i ) 


r 2 

r i+l/2 





(CI) 



r 



, (^+3/2 ^+1/2) 4 (p'i+1/2 ^1-1/2 J 



i+i ' 



( r i+3/2 ^+1/2) {'"I 



(+1/2 r f-l/2. 



while to evolve B we difference equation (^) according to 



?'< 1 _ J?" ^f^l'/T _ - -" ' 



B'.' V '+1/2 s i-l/2; 

= rA, — -r- (C2) 



A/ '* ^.2 



W/2~ r *-l/2 



Boundary condition ( |l2[ ) at the origin requires 



Bf 1 = 0. (C3) 



The outer boundary condition depends on the nature of the exterior medium: 

B"^ x = 0, (vacuum exterior) (C4) 

according to equation (|l3]), or 

B" +l =B" -. — \\T^](B" A, (plasma exterior), (C5) 



according to equation (|l4|). The two prescriptions agree in the case of a vacuum exterior, where p ex = and 1Z = -1, whereby 
condition (]C5|) gives Bf""= fi" =,..., = B? =0 



Substituting finite difference equation ( p2| ) into (CI) and using equation (C3) at the inner boundary and equation (24) or (25) at 
the outer boundar y yi elds a single tridiagonal matrix equation for fl n+l . We invert this tridiagonal matrix at each timestep, after which 
we use equation ( |C2[ ) again to get B n+1 . Since we use implicit differencing in time, the upper bound on the size of our integration 
timestep At is set by accuracy and not stability criteria. 
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